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I.  INTRODUCTION 

A.      MISSION  NEED 

Unmanned  Aerial  Vehicles  (UAVs)  are  gaining  acceptance  as  an  integral  part  of  the 
operations  of  today's  armed  forces.  Preceding  and  during  Operation  Desert  Storm,  UAVs 
flew  a  variety  of  missions  including  reconnaissance,  targeting  for  gunfire  support,  and 
battle  damage  assessment.  Desert  Storm  provided  the  first-ever  combat  test  of  Unmanned 
Aerial  Vehicles  by  U.S.  forces  [Ref.  I].  Reviews  of  lessons  learned  from  Desert  Shield 
and  Desert  Storm  reveal  the  outstanding  successes  of  UAVs. 

There  are  many  examples  of  effective  UAV  use  during  the  war.  In  one  instance,  a 
commander  of  a  Marine  task  force  was  able  to  monitor  UAV  imagery  of  Kuwait  as  his 
task  force  approached  the  city,  revealing  the  exact  reaction  of  the  Iraqi  forces  to  Marine 
armor,  artillery,  and  troop  movements.  The  Navy  used  UAVs  to  search  for  mines,  spot 
for  gunfire  support,  perform  reconnaissance  missions  for  SEAL  teams,  and  search  for 
Iraqi  Silkworm  sites,  command  and  control  bunkers,  and  anti-aircraft  artillery  sites.  The 
Marines  quickly  reacted  to  an  Iraqi  attack  into  Saudi  Arabia  observed  by  a  UAV  and 
decisively  crushed  the  Iraqi  invasion  with  airborne  Cobras  and  Harriers    The  Army 
provided  their  Apache  pilots  with  route  reconnaissance  acquired  from  UAVs  shortly 
before  the  Apache  missions.  Spotting  for  air  strikes  and  naval  gunfire  support 


became  so  successful  that  Iraqi  soldiers  were  seen  attempting  to  surrender  to  UAVs  as 
they  flew  overhead  [Ref.  2]. 

UAVs  have  many  advantages  over  manned  aircraft  which  help  account  for  their 
effectiveness.  First,  the  cost  of  a  UAV  is  a  very  small  fraction  of  the  cost  of  a  manned 
aircraft.  The  Pioneer,  for  example,  costs  approximately  $500,000  for  the  aircraft  and 
$500,000  for  the  onboard  camera,  bringing  the  total  cost  of  the  package  to  a  mere  one  to 
two  percent  of  the  cost  of  most  manned  tactical  aircraft.  UAVs  are  also  extremely  flexible 
with  respect  to  launch  platform.  They  can  be  launched  from  small  fields,  truck  beds,  and 
practically  any  ship  in  the  Navy  inventory.  UAVs  are  frequently  very  hard  to  detect  with 
radar  or  infrared  systems  due  to  their  small  size,  composite  construction,  small  engines 
(sometimes  electric  motors),  and  their  slow  speed.  UAVs  also  have  an  advantage  from 
being  unmanned.  The  aircraft  is  not  limited  by  the  "g"  tolerance  of  a  pilot  or  by  pilot 
fatigue.  Finally,  the  best  advantage  of  all  is  that  when  a  UAV  crashes  or  is  lost  to  enemy 
fire,  there  is  no  search  and  rescue  mission  required,  no  prisoners  of  war  taken,  and  no  loss 
of  life. 

UAVs  are  not  without  their  problems,  however.  Acquisition  and  support  programs 
are  relatively  new  and  underdeveloped.  This  results  in  a  UAV  force  that  is  relatively  small 
in  number  of  aircraft,  and  small  and  inexperienced  in  terms  of  personnel.  The  Pioneer 
showed  signs  of  these  underlying  problems  during  Desert  Storm.  Six  Pioneers  were 
damaged  badly  enough  to  require  return  to  the  factory  for  repair  due  to  no  intermediate 
level  maintenance  facilities  being  available.  Five  Pioneers  were  lost  due  to  mechanical 


malfunction  or  operator  error  [Ref.  2].  Due  to  the  small  number  of  available  aircraft, 
losses  are  very  costly  and  operator  errors  need  to  be  avoided  if  possible.  Thorough  and 
frequent  training  through  simulation  can  keep  operators  performing  at  their  peak. 
B.       STATEMENT  OF  OBJECTIVE 

To  provide  realistic  simulation  for  training,  accurate  aerodynamic  data  for  the 
aircraft  must  be  available.  Aerodynamic  parameters  of  the  aircraft  can  be  determined  from 
its  physical  characteristics,  wind  tunnel  tests,  and  flight  tests  of  actual  or  scale  models.  It 
is  the  objective  of  this  work  to  provide  the  ground  work  for  determining  whether  flight 
test  can  effectively  be  used  to  accurately  estimate  the  static  and  dynamic  stability  and 
control  characteristics  of  a  fixed  wing  UAV.  Flying-qualities  parameters  were  estimated 
using  analytic  techniques,  and  the  resultant  flight  dynamics  were  simulated  to  provide 
expected  behavior  for  future  test  flights. 


H.  BACKGROUND 

A.      EVOLUTION  OF  AN  AIRCRAFT 

For  any  aircraft  to  get  from  an  idea  to  an  actual  flying  vehicle,  properties  such  as 
stability  and  control  derivatives,  which  determine  the  flying  characteristics  of  the  aircraft, 
must  be  determined.  These  derivatives  are  used  to  size  control  surfaces,  design  flight 
control  systems,  and  program  training  devices  such  as  simulators.  There  are  typically 
three  ways  to  determine  or  estimate  derivatives. 

The  first  method  of  determining  an  aircraft's  derivatives  begins  early  and  continues 
throughout  the  design  process.  This  step  involves  determining  derivatives  mathematically 
from  physical  characteristics  of  the  aircraft.  Requiring  little  more  than  a  few 
well-educated  engineers  and  some  calculators  or  desktop  computers,  this  method  is 
relatively  simple.  Derivatives  can  be  reasonably  approximated,  but  must  be  refined 
through  other  methods. 

The  second  method  of  determining  derivatives  is  through  aerodynamic 
measurements  of  wind  tunnel  testing.  This  method  is  more  complicated  than  simple  pencil 
and  paper  calculations  due  to  the  necessity  for  model  making  and  wind  tunnel  operation, 
but  much  better  results  can  be  achieved.  Derivatives  must  still  be  refined,  however,  due  to 
factors  such  as  scale  effects  and  interference  from  wind  tunnel  walls  and 


supporting  hardware.  Also,  dynamic  effects  are  often  difficult  to  properly  account  for  in 
most  wind  tunnels. 

The  final  step  approach  to  determining  the  derivatives  of  an  aircraft  is  through  flight 
test.  This  method  is  by  far  the  most  expensive  and  complicated  way  to  determine  the 
aircraft's  derivatives,  but  it  is  also  the  most  precise  and  complete.  Scaled  flight  test 
provides  a  viable  option  for  this  third  method. 
B.       PIONEER  UAV 

In  June  1982,  Israeli  forces  very  successfully  used  UAVs  as  a  key  element  in  their 
attack  on  Syria.  Scout  and  Mastiff  UAVs  were  used  to  locate  and  classify  SAM  and  AAA 
weaponry  and  to  act  as  decoys  for  other  aircraft.  This  action  resulted  in  heavy  Syrian 
losses  and  minimal  Israeli  losses.  A  year  and  a  half  later,  the  U.S.  Navy  launched  strikes 
against  Syrian  forces  in  the  same  area  with  losses  much  higher  for  the  Navy  than  those  of 
the  Israelis  [Ref.  3]. 

The  Commandant  of  the  Marine  Corps,  General  P.  X.  Kelly,  recognized  the 
effectiveness  of  the  Israeli  UAVs.  Secretary  of  the  Navy  John  Lehman  then  initiated 
development  of  a  UAV  program  for  the  U.S.  Anxious  to  get  UAVs  to  the  fleet,  Secretary 
Lehman  stipulated  that  UAV  technology  would  be  off-the-shelf  [Ref.  4].  After  the 
contract  award  to  AAI  Corporation  of  Baltimore,  Maryland  for  the  Pioneer  UAV, 
developmental  and  operational  testing  took  place  concurrently.  This  approach  resulted  in 
quick  integration  of  the  Pioneer  into  the  fleet.  Unfortunately,  such  quick  integration  into 


the  fleet  can  result  in  problems  identified  during  operational  use  which  had  not  been  fully 
explored  in  the  test  and  evaluation  process. 

The  UAV  Office  at  the  Pacific  Missile  Test  Center  (PMTC,  now  the  Naval  Air 
Warfare  Center,  NAWC,  Weapons  Division,  Pt.  Mugu)  was  tasked  with  Developmental 
Test  and  Evaluation  of  the  Pioneer.  Testing  revealed  the  following  concerns  which 

warranted  further  investigation  [Refs.  5,6]: 

1 .  discrepancies  in  predicted  with  flight-tested  rate  of  climb,  time  to  climb, 
and  fuel  flow  at  altitude; 

2.  apparent  autopilot-related  pitch  instability; 

3 .  tail  boom  structural  failure; 

4.  severely  limited  lateral  control; 

5.  slow  pitch  response  causing  degraded  maneuverability  at  high  gross 
weights; 

6.  insufficient  testing  to  determine  the  effects  of  the  new  wing  on  flight 
endurance. 

The  Target  Simulation  Laboratory  at  Pt.  Mugu  was  tasked  to  develop  a  computer 
simulation  of  the  Pioneer  in  order  to  provide  cost-effective  training  for  pilots. 
Aerodynamic  data  were  needed  to  provide  the  stability  and  control  derivatives  necessary 
for  the  simulation  as  well  as  to  answer  questions  concerning  basic  flying  qualities  of  the 
Pioneer. 

In  order  to  provide  support  to  the  research  being  done  at  Pt.  Mugu  and  to  provide 
for  future  UAV  project  support,  a  research  program  was  begun  at  the  Naval  Postgraduate 
School  (NPS).  An  instrumented  half-scale  radio-controlled  model  of  the  Pioneer  was  used 


for  the  research  at  NPS.  Research  performed  included  wind  tunnel  tests,  flight  tests,  and 
numerical  modeling. 

Initial  NPS  research  on  the  Pioneer,  performed  by  Capt.  Daniel  Lyons,  involved  a 
computer  analysis  of  the  Pioneer  in  its  original  configuration  and  with  a  proposed  larger 
tail.  A  low  order  panel  method  (PMARC)  was  used  for  the  aerodynamic  analysis.  Static 
longitudinal  and  directional  stability  derivatives,  the  neutral  point,  and  crosswind 
limitations  were  calculated.  Drag  polars  were  constructed  using  the  component  buildup 
method  for  profile  drag,  and  drag  reduction  measures  were  considered  [Ref  7]. 

In  conjunction  with  Capt.  Lyons  work,  Lt.  James  Tanner  conducted  wind  tunnel 
tests  to  determine  propeller  efficiencies  and  thrust  coefficients  for  drag  studies  [Ref.  8]. 
Lt.  Tanner  also  conducted  flight  tests  to  determine  power  required  curves  and  drag  polars 
[Ref.  8].  Capt.  Robert  Bray  later  conducted  wind  tunnel  tests  of  a  0.4-scale  model  at 
Wichita  State  University  to  determine  static  stability  and  control  derivatives  [Ref.  9]. 
Aerodynamic  data  obtained  by  Capt.  Lyons  and  Bray  have  been  supplied  to  PMTC  to  be 
used  for  simulation. 

Lt.  Jim  Salmons  performed  initial  flying  qualities  flight  testing  using  an  onboard  data 
recording  system  in  order  to  determine  static  stability  parameters.  Unfortunately, 
vibration  problems  with  the  onboard  recorder  rendered  much  of  the  data  unusable  [Ref. 
10]. 

Following  up  on  Lt.  Salmons'  work,  Lt.  Kent  Aitcheson  installed  the  CHOW-1G 
telemetry  system,  designed  by  Lt.  Kevin  Wilhelm,  in  an  attempt  to  alleviate  the  vibration 


problem  experienced  by  Lt.  Salmons.  The  new  flight  test  configuration  was  used  to  test 
static  longitudinal  and  lateral-directional  stability  characteristics  of  the  Pioneer.  The 
vibration  problem  experienced  by  Lt.  Salmons  was  overcome,  though  not  enough  data 
were  acquired  for  a  complete  and  thorough  analysis  of  the  Pioneer's  characteristics.  Much 
insight  was  gained,  however,  concerning  instrumentation.  Resolution  needed  to  be 
improved  for  flight  control  position  indication  [Refs.  11,12]. 

Lt.  Paul  Koch  conducted  further  flight  tests  of  the  Pioneer  with  the  CHOW-1G 
telemetry  system.  Static  longitudinal  stability  results  from  the  flight  tests  correlated  well 
with  theoretical  predictions  and  with  simulations  of  a  full-scale  Pioneer.  Electromagnetic 
interference  with  the  flight  control  system  at  the  test  site  resulted  in  loss  of  the  half-scale 
model  Pioneer  before  further  data  collection  and  analysis  could  be  performed  [Ref  13]. 
C.       PARAMETER  ESTIMATION 

Parameter  estimation  is  used  to  derive  stability  and  control  derivatives  from  dynamic 
flight  test  data.  Lcdr.  Robert  Graham  successfully  used  the  Modified  Maximum 
Likelihood  Estimation  (MMLE)  technique  with  MATLAB  software  to  analyze  simulated 
and  actual  flight  test  data  of  several  aircraft.  Analysis  of  simulated  UAV  data  revealed  the 
effect  of  signal-to-noise  ratio  on  the  estimator,  and  the  need  for  proper  control-surface 
excitation  for  a  particular  response  [Ref.  14].  Cdr.  Patrick  Quinn  analyzed  flight  test  data 
from  the  Marine  Corps  BQM-147  UAV  using  both  MMLE  and  a  more  robust  non-linear 
model,  pEst,  and  compared  the  results  of  the  two  approaches.  Noise  was  found  to  be  a 
problem  when  the  system  response  was  in  the  same  general  frequency  as  the  noise. 


Limited  available  data  and  noise  at  the  frequency  expected  for  the  system's  response 
prevented  a  successful  resolution  of  all  stability  and  control  derivatives  of  interest.  The 
pEst  model  was  found  to  be  the  estimator  of  choice  when  aircraft  maneuvers  exceed  what 
is  generally  considered  reasonable  for  linear  approximation  of  flight  dynamics  [Ref.  15]. 
While  previous  work  with  UAVs  at  NPS  has  at  times  been  frustrating,  much  has 
been  learned,  especially  concerning  the  most  challenging  method  of  aircraft  analysis,  flight 
test.  This  work  initiates  the  implementation  of  the  lessons  learned  from  previous  work 
into  dynamic  simulation  and  flight  testing  of  a  generic  UAV  in  order  to  properly  prepare 
the  UAV  lab  at  NPS  for  further  support  of  future  projects  and  to  demonstrate  the  value  of 
scaled  UAV  flight  test  to  the  fleet. 


m.  THE  AIRCRAFT 

A.  GENERAL  DESCRIPTION 

The  "Bluebird",  shown  in  Figures  1-3,  is  a  high-wing  tricycle-gear  radio-controlled 
airplane.  It  is  constructed  of  wood,  foam,  composites,  and  metal.  It  is  powered  by  a 
Sachs-Dolmar  3.7-cubic-inch  two-stroke  gasoline  engine  which  drives  a  24  inch 
two-bladed  wood  propeller.  It  is  controlled  by  a  nine-channel  pulse-code-modulated 
Futaba  radio  operating  at  72.710  MHz.  To  enhance  reliability,  the  Bluebird  has  two 
receivers  which  share  control  of  the  aircraft.  The  left  receiver  controls  the  left  aileron, 
elevator,  and  flap,  and  engine  ignition  and  onboard  electronics  package  cut-off,  while  the 
right  receiver  controls  the  right  aileron,  elevator,  and  flap,  and  rudder,  nose-wheel 
steering,  and  throttle.  The  Bluebird  can  fly  within  visual  range  for  approximately  1.5 
hours.  Table  1  describes  physical  specifications  of  the  Bluebird. 

B.  STABILITY  DERIVATIVES 

The  initial  estimates  of  the  stability  derivatives  of  the  Bluebird  were  made  using  the 
physical  characteristics  of  the  aircraft  such  as  airfoil  data  ,  geometric  measurements, 
relative  positions  of  aircraft  components,  mass,  and  weight  [Refs.  16-19].  The  assumed 
flight  condition  with  the  associated  aircraft  configuration  is  described  in  Table  2. 
Nondimensional  stability  and  control  derivatives  estimated  are  shown  in  Table  3,  and 
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dimensional  stability  and  control  derivatives  estimated  are  shown  in  Table  4.  MATLAB 

programs  written  for  the  physical  and  derivative  calculations  appear  in  Appendix  A. 

TABLE  1 
SPECIFICATIONS 

Length   9.84  ft. 

Height    3.04  ft. 

Wing  Airfoil  (est.)    GO  769 

Horizontal  Stab.  Airfoil  (est.)    NACA4412 

S.ng(Sref)    22.38ft. 

S 4.701  ft.2 


2 


Sv  1.277  ft.2 

c  1.802  ft. 

ct    1.281  ft. 

b  12.42  ft. 

b,    3.67  ft. 

bv   1.21ft. 

AR  6.89 

AR,    2.86 

AR,    1.14 

VH    0.6274 

Vv 0.0247 

V   0.0023 

lt   5.381  ft. 

Iv  5.381  ft. 


TABLE  2 
FLIGHT  CONDITION/AIRCRAFT  CONFIGURATION 

Weight  57.79  lbs. 

!„   12.58  slug*ft2 

lyy   13.21  slug*ft2 

I,,   19.99  slug*ft2 

Velocity  60  mph/88  ft/sec 

Altitude    800  ft  MSL 

Density    0.002327  slugs/ft3 

Center  of  Gravity  27%  of  m.a.c. 

CUrim    0.2866 

A.o.A,rim   3.8  degrees 

Elevatorlrjm    1.6  degrees 
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Figure  1 
Side  View 


Figure  2 
Top  View 
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Figure  3 
Front  View 
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TABLE  3 
NONDIMENSIONAL  DERIVATIVES 


CL    0.2866 

Cu    4.1417 

CL.    1.5787 

C„«    -1.0636 

CLq    3.9173 

CLhe   0.4130 

Cmu    -1.2242 

CH    0.0484 

Cyp     0.0000 

Ci,   -0.3579 

Cn -0.0526 

Cysa   0.0000 

Ciu    0.2652 

Cnhr    -0.0326 

CDq    0.0000 


CD    0.0358 

CDo    0.0311 

CDa    0.1370 

Cm.    -4.6790 

Cm]    -11.6918 

CDu  0.0650 

Cy9    -0.3100 

C/p    -0.0330 

Cnp     -0.0358 

Cyr   0.0967 

Ci 0.0755 

Cnsa    -0.0258 

Cyt,    0.0697 

C/5r     0.0028 

Cyq   0.0000 


TABLE  4 
DIMENSIONAL  DERIVATIVES 


Xu    -0.0914 /sec 

Xse    -7.2961  ft/sec2 

Za    -468.9852  ft/sec2 

Zq    4.5027  ft/sec 

Mu  0.0000  /ft*sec 

M-a    -1.3178 /sec 

Mse    -33.6730  /sec2 

Yp    0.0000  ft/sec 

Ysa  0.0000  ft/sec2 

Lp    -6.5787  /sec2 

Lr    1.0613 /sec 

Lsr  0.5589  /sec2 

Np    -0.3167 /sec 

Nsa    -3.2375  /sec2 


Xa    16.7894  ft/sec 

Zu  -0.7312  /sec 

Z-a    1.8146  ft/sec 

lie   -46.368  ft/sec2 

Ma  -29.2559  /sec2 

Mq   -3.2928  /sec 

7p  -34.8021  ft/sec2 

Yr  0.7663  ft/sec 

YSr    -7.8282  ft/sec2 

Lp  -5.0281  /sec 

Lsa   52.7966  /sec2 

A^p    6.0593  /sec2 

Nr    -0.4647  /sec 

A^sr  -4.0900  /sec2 
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C.      MOMENTS  OF  INERTIA 

Having  accurate  moments  of  inertia  is  critical  to  ensuring  the  accurate  prediction  of 
aircraft  dynamics.  Direct  calculation  of  a  model's  moments  of  inertia  by  consideration  of 
the  contributions  made  by  individual  parts  is  impractical  and  inaccurate.  Determination  of 
the  moments  of  inertia  by  test  is  much  more  practical  and  precise.  The  changes  in  a 
model's  moments  of  inertia  due  to  addition  or  subtraction  of  equipment  or  structure  can  be 
calculated  directly,  thereafter. 

In  the  determination  of  moments  of  inertia  by  test,  the  aircraft  is  hung  from  the 
ceiling  and  swung.  Using  the  period  exhibited  and  the  principles  of  compound  pendulums, 
the  moments  of  inertia  of  the  model  can  be  extracted  [Ref  20-21].  In  order  to  calculate 
the  moment  of  inertia  about  all  axes,  the  model  must  be  hung  from  the  ceiling  and  swung 
three  different  ways,  each  such  that  as  the  aircraft  swings,  it  is  rotating  about  the  axis  of 
interest.  The  Bluebird  was  hung  by  chain  and  swung  as  pictured  in  Figures  4-6. 
Specifications  for  the  geometry  of  each  test  can  be  found  in  Appendix  B. 

Reference  2 1  provides  equation  ( 1 )  for  calculating  the  moment  of  inertia  of  a 

swinging  model. 

j  -  *L 

4n2    fm+s        4n2ys         g 


j         Wm+sZm+s  Jl  WsZs  Jl         WMZ2M  nv 

7  =       ..2    Pm+s  ~  -TTPs  ~  -g—  V> 


W  is  the  weight,  Z  is  the  distance  from  pivot  to  center  of  gravity,  p  is  the  period,  and  g  is 
the  gravitational  constant.  M  and  S  are  subscripts  designating  either  the  model  or  the 

support. 
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Figure  4 
1^  Test  (not  to  scale) 
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Figure  5 
Iyy  Test  (not  to  scale) 
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Figure  6 
Ia  Test  (not  to  scale) 
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It  was  determined  that  swinging  the  support  (the  chains)  in  the  configuration  it 
would  be  in  when  supporting  the  model  would  not  be  possible,  since  the  chains  would  not 
maintain  their  positions  without  the  model  in  place.  Equation  (1)  was  therefore 
manipulated  in  order  to  treat  the  chains  as  long  slender  rods  and  to  calculate  their 
moments  of  inertia  as  such  [Ref.  22]. 

The  new  form  of  equation  (1)  is 

7  =       4,1    Pm+s  ~  -g~  ~  S      3g  (2) 

where  Ls  is  the  length  of  a  chain  and  the  summation  is  taken  over  all  chains  (four  in  this 
case).  In  particular,  there  were  two  "long"  chains  and  two  "short"  chains,  all  having  a 
weight  per  unit  length  co .  Appropriate  substitutions  were  made  to  yield 

T          Wm  +  2a(LsHORT  +  LlqngWm+s  n                 iVMZu  2a(T2  .    Tl         \  (X\ 

1    -    ^3 ^M+S g -^\L  SHORT +L  LONG  J  V    ' 

Having  the  equation  in  this  form  fixes  the  values  for  all  variables  except 
Z\i+s,  Zm,  and  Pm+s.  These  three  variables  were  measured  for  each  of  the  three 
configurations  and  calculations  were  made.  Four  periods  were  timed  during  the  swing 
tests,  and  the  tests  were  repeated  ten  times.  The  moments  of  inertia  calculated  are  shown 
in  Table  2. 
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IV.  PREDICTED  DYNAMICS 

A.  GENERAL 

When  flight  testing  an  aircraft,  it  is  important  to  attempt  to  predict  the  results  of  the 
testing  prior  to  the  actual  flights.  The  information  provided  by  the  predictions  can  be  used 
in  briefing  the  pilot  as  to  what  to  expect  from  the  aircraft  and  also  to  avoid  any  potentially 
dangerous  flight  regimes.  Longitudinal  and  lateral-directional  dynamics  of  the  Bluebird 
have  been  predicted  using  the  stability  and  control  derivatives  estimated  in  chapter  three 
along  with  computational  methods  based  upon  the  six  equations  of  motion  as  described  in 
references  18  and  23.  Computational  programs  were  written  in  MATLAB  and  appear  in 
Appendix  C. 

B.  LONGITUDINAL  DYNAMICS 

The  MATLAB  program  named  "longnat.m"  uses  the  full  (4x4)  longitudinal  plant 
and  the  MATLAB  "impulse"  function  to  determine  the  short  and  long  period  natural 
responses.  A  diary  of  the  values  computed  and  output  by  the  program  provides 
eigenvalues,  damping  ratios,  and  damped  and  undamped  natural  frequencies  for  both 
modes  as  shown  in  Table  5. 
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TABLE  5 
LONGITUDINAL  DATA 

Short  Period  Long  Period 

Eigenvalues                                     -5.083  +/-  4.861  i  -0.037  +/-  0.400i 

Damping  ratio                                       0.723  0.093 

Undamped  Natural  Frequency         7.03  rad/sec  0.401  rad/sec 

Damped  Natural  Frequency            4.86  rad/sec  0.399  rad/sec 

Figures  7  and  8  show  the  combined  short  and  long  period  natural  response.  It  can  be  seen 
that  the  short  period  response  is  almost  completely  damped  out  after  only  two  seconds. 
Response  beyond  two  seconds  is  primarily  long  period.  The  phugoid  mode  is  seen  to  be 
lightly  damped. 

The  short  and  long  period  response  to  a  unit  step  elevator  input  is  determined  by 
using  the  full  (4x4)  longitudinal  plant  and  the  MATLAB  "step"  function  in  the  MATLAB 
program  named  "  stepper,  m".  Figures  9  and  10  show  the  short  and  long  period  response 
to  a  step  input.  In  the  first  plot,  the  long  period  can  be  seen  to  be  much  more  heavily 
excited  than  the  short  period  mode.  Again,  it  can  be  seen  that  the  short  period  response  is 
almost  completely  damped  out  after  about  two  seconds.  Due  to  held  elevator  input,  the 
long-term  response  is  to  trim  to  a  new  angle  of  attack,  while  pitch  rate  dies  out. 

The  MATLAB  program  named  "n_step.m"  uses  the  full  (4x4)  longitudinal  plant  and 
the  MATLAB  "step"  function  to  determine  the  normal  acceleration  or  load  factor 
response  to  a  unit  step  elevator  input.  Figure  1 1  shows  the  short  period  normal 
acceleration  response  to  a  unit  step  input.  It  should  be  noted  that  normal  acceleration  is 
defined  opposite  in  sign  to  load  factor.  The  long-period  response  is  seen  to  be  much 
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Figure  7 
Natural  Response  (short  time) 


Figure  8 
Natural  Response  (long  time) 
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Response  to  Step  Input  (short  time) 


0 


O 


0) 

2-1 
o 

a. 


-3 1 


-5 


_...a 

r                                      1 

T\ 

<      i 

—j \ 

!   1       \      t 

ir             t         : 

fx k 

^^"^~<- — 

i  ■        '  / 

"alpha  ~X~ 

~i 

1 1 
1 1 

I* 

V 

i 

i 

20 


40  60 

Time,  sec 


80 


100 


Figure  10 
Response  to  Step  Input  (long  time) 
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more  critical  than  the  short-period  response  in  generating  large  load  factors  for  this 
aircraft. 


Figure  1 1 
Normal  Acceleration  Response  to  Step  Input 

The  longitudinal  response  to  given  initial  conditions  is  determined  using  the  full 
(4x4)  longitudinal  plant  and  the  MATLAB  "initial"  function  in  the  MATLAB  program 
named  "homogen.m".  The  initial  conditions  (  u/U  =  .34,  alpha  =  5  deg,  q  =  8.8  deg/sec, 
theta  =  -0.8  deg)  are  provided  for  all  states  from  the  step  input  results  after  initial 
dynamics  have  died  out  (approximately  15  seconds).  Figures  12  and  13  show  the 
longitudinal  response  to  the  initial  conditions.  The  response  in  angle  of  attack  is  small 
while  a  significant  pitch  rate  is  developed  in  both  the  short-period  and  long-period 
responses.  The  long-period  is  again  seen  to  be  lightly  damped. 

The  MATLAB  program  named  "doublet. m"  uses  the  full  (4x4)  longitudinal  plant 
and  the  transition  matrix  (using  the  MATLAB  "exponential"  function)  to  find  the 
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Figure  12 
Homogeneous  Response  (short  time) 
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Homogeneous  Response  (long  time) 
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response  to  a  unit  pitch  doublet  input  at  approximately  the  short  period  damped  natural 
frequency.  Figure  14  shows  the  input  pitch  doublet,  and  Figure  15  shows  the  response. 
Though  the  doublet  is  commonly  used  to  observe  only  the  short-period  response,  the 
angle  of  attack  response  indicates  the  long-period  response  is  also  excited. 

The  response  (transfer  function  gain  and  phase)  to  a  harmonic  elevator  input  is 
found  using  the  full  (4x4)  longitudinal  plant  and  the  MATLAB  "bode"  function  in  the 
MATLAB  program  named  "sp_bode.m".  Figures  16  and  17  show  the  gains  and  phase 
shifts  from  an  elevator  frequency  sweep.  The  angle  of  attack  gain  can  be  seen  decreasing 
from  a  peak  at  a  very  low  excitation  frequency;  the  short-period  response  is  masked  by  the 
long-period  gain.  In  actuality,  the  low-frequency  response  is  of  little  interest.  Pitch  rate 
shows  a  maximum  gain  at  approximately  2.4  radians/second  corresponding  to  an  elevator 
cycle  period  of  approximately  2.6  seconds. 

The  MATLAB  program  named  "nbode.m"  uses  the  full  (4x4)  longitudinal  plant 
and  the  MATLAB  "bode"  function  to  find  the  normal  acceleration  or  load  factor  response 
(transfer  function  gain  and  phase)  to  a  harmonic  input.  Figures  18  and  19  show  the 
normal  acceleration  gains  and  phase  shifts  from  an  elevator  frequency  sweep.  The  short 
and  long-period  damped  natural  frequencies  can  be  observed  at  4.86  radians/second  and 
0.399  radians/second  respectively  where  their  respective  gains  peak.  While  the  gain 
demonstrated  at  the  long-period  damped  natural  frequency  is  quite  significant,  the  aircraft 
would  not  be  expected  to  be  excited  at  that  frequency.  An  excitation  near  the 
short-period  damped  natural  frequency  is  much  more  likely.  Flight  test  pilots  and 
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Figure  14 
Unit  Elevator  Doublet  Input 
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Figure  16 
Frequency  Response  (Gain) 
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Figure  17 
Frequency  Response  (Phase) 
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Figure  18 
Normal  Acceleration  Response  to  a  Harmonic  Input  (Gain) 


Figure  19 
Normal  Acceleration  Response  to  a  Harmonic  Input  (Phase) 
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engineers  should  be  aware  of  the  fact  that  excitation  at  that  frequency  will  produce 

approximately  0. 15  "g's"  per  degree  of  elevator  deflection  or  that  it  takes  a  harmonic 

amplitude  of  about  6.7  degrees  of  elevator  to  produce  one  "g"  of  acceleration.  The 

Bluebird  has  a  maximum  of  15  degrees  of  up  elevator  and  12  degrees  of  down  elevator 

available.  This  would  equate  to  approximately  -0.8  to  +3.2  "g's",  which  is  easily  tolerated 

structurally. 

C.      LATERAL-DIRECTIONAL  DYNAMICS 

Undamped  natural  frequency,  damped  natural  frequency,  damped  natural  period, 
and  damping  ratio  for  the  Dutch  roll  mode;  time  constant  for  the  roll  mode;  and  time 
constant  and  time  to  double  or  half  for  the  spiral  mode  are  calculated  by  the  MATLAB 
program  named  "latdir.m"  using  the  full  (4x4)  lateral-directional  plant.  Values  calculated 

are  described  as  follows. 

Dutch  roll  mode: 

Undamped  natural  frequency  2.65  rad/sec 

Damped  natural  frequency  2.62  rad/sec 

Damped  natural  period  2.40  sec 

Damping  ratio  0.148 

Roll  mode: 

Time  constant  0.195  sec 

Spiral  mode: 

Time  constant  -29.28  sec 

Time  to  double  amplitude  20.29  sec 

The  Dutch  roll  mode  is  observed  to  be  lightly  damped.  Also,  the  spiral  mode  is  divergent 
with  a  fairly  short  time  to  double  bank  angle. 
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The  MATLAB  program  named  "rudkick.m"  uses  the  full  (4x4)  lateral-directional 
plant  and  the  MATLAB  "impulse"  function  to  find  the  response  to  a  unit  rudder  impulse. 
The  program  also  finds  the  bank  angle  at  the  end  of  100  seconds  which  is  approximately 
four  degrees  in  the  direction  of  the  rudder  kick.  Figure  20  shows  the  unit  rudder  kick 
response  in  bank  and  sideslip  angles.  Sideslip  lags  bank  angle  by  approximately  0.6 
seconds  or  80  degrees  of  phase  and  is  approximately  60  percent  larger  in  magnitude. 

The  response  to  an  initial  sideslip  is  determined  by  the  MATLAB  program  named 
"sideslip. m"  using  the  full  (4x4)  lateral-directional  plant  and  the  MATLAB  "initial" 
function.  Initial  conditions  for  sideslip  and  bank  angle  are  provided  by  a  steady-sideslip 
condition  where  bank  angle  is  ten  degrees,  and  sideslip  angle  is  13.7  degrees.  Figure  21 
shows  the  homogeneous  response  to  the  initial  sideslip.  The  lightly  damped  Dutch  roll 
mode  can  be  seen  by  observing  the  oscillations  of  both  bank  angle  and  sideslip  angle,  while 
the  divergent  spiral  mode  can  be  observed  by  the  increasing  bank  angle. 

The  MATLAB  program  named  "roll.m"  uses  the  full  (4x4)  lateral-directional  plant 
and  the  MATLAB  "bode"  function  to  find  the  roll  angle  response  (transfer  function  gain 
and  phase)  due  to  a  harmonic  aileron  input.  Figures  22  and  23  show  the  roll  angle  gains 
and  phase  shifts  from  the  harmonic  aileron  input.  The  maximum  high-frequency  roll  angle 
gain  is  seen  to  occur  at  approximately  2.8  radians/second  which  corresponds  to  an  aileron 
input  period  of  approximately  2.2  seconds. 
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Figure  20 
Unit  Rudder  Kick  Response 
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Figure  21 
Homogeneous  Response  to  an  Initial  Sideslip 
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Figure  22 
Roll  Angle  Gain  due  to  a  Harmonic  Aileron  Input 


Figure  23 
Roll  Angle  Phase  Shift  due  to  a  Harmonic  Aileron  Input 
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D.      REMARKS 

Analysis  has  been  conducted  for  the  most  common  modes.  The  Bluebird 
demonstrates  no  particularly  dangerous  flying  qualities.  Particular  characteristics  of  the 
Bluebird's  behavior  worth  noting  include  a  very  heavily  damped  short  period  longitudinal 
mode,  a  lightly  damped  phugoid  mode,  a  lightly  damped  Dutch  roll  mode,  and  a  divergent 
spiral  mode.  Additional  analysis  could  be  easily  done  by  analogy  to  those  modes  all  ready 
analyzed.  Further  analysis  might  include,  for  example,  the  response  to  a  step  aileron 
input,  aileron  doublet,  or  elevator  impulse  ("stick  rap"). 
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V.  FLIGHT  TEST 

A.       TEST  DESCRIPTION 

The  relative  positions  of  an  aircraft's  center  of  gravity  and  neutral  point  are  critical 
to  the  aircraft's  longitudinal  stability  and  handling  qualities.  In  order  for  a  conventional 
(aft  tail)  aircraft  to  be  statically  stable,  its  center  of  gravity  must  be  forward  of  its  neutral 
point.  The  farther  the  center  of  gravity  is  in  front  of  the  neutral  point  (relative  to  the 
chord  length),  the  more  statically  stable  the  aircraft  is.  As  the  center  of  gravity  is  moved 
aft  beyond  the  neutral  point,  the  aircraft  becomes  statically  unstable  and  more 
maneuverable. 

As  an  aircraft's  center  of  gravity  is  moved  aft  toward  the  neutral  point,  less  and  less 
change  in  elevator  trim  is  required  to  achieve  steady,  level  flight  for  a  given  change  in 
airspeed.  This  fact  can  be  used  to  experimentally  determine  the  aircraft's  neutral  point. 
Reference  23  provides  equation  (4)  which  describes  the  relationship  between  change  in 
pitching  moment  coefficient  with  change  in  coefficient  of  lift  and  the  required  change  in 
elevator  deflection  with  change  in  coefficient  of  lift. 

dCL  H     ^L'i*      dCL  v   ' 

As  the  center  of  gravity  moves  aft  and  approaches  the  neutral  point,  dCm/dCL  approaches 
zero.  Since  VH  and  CL5e  are  constants,  d6e/dCL  must  also  approach  zero. 
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In  order  to  determine  the  neutral  point  by  flight  test,  the  aircraft  must  be  flown  at 
various  centers  of  gravity.  At  each  particular  center  of  gravity,  the  aircraft  is  trimmed  at 
several  airspeeds  (coefficients  of  lift).  Each  airspeed  and  the  corresponding  elevator 
deflection  are  recorded.  Plots  are  then  generated  showing  elevator  deflection  versus 
coefficient  of  lift  at  each  center  of  gravity  tested.  A  line  is  then  best  fit  through  the  data 
points  for  each  center  of  gravity.  The  slope  of  this  line  represents  d5e/dCL  for  that 
particular  center  of  gravity.  Finally,  d5e/dCL  versus  center  of  gravity  is  plotted.  A  best  fit 
line  is  then  drawn  through  these  data  points.  Using  the  best  fit  line  just  drawn,  the  center 
of  gravity  where  d5e/dCL  equals  zero  is  the  aircraft's  neutral  point. 
B.       INSTRUMENTATION 

In  order  to  acquire  the  airspeed  of  the  Bluebird,  a  simple,  commercially-available 
airspeed  indicator  was  installed.  The  Digicon  TT-01  Tele  Tachometer/ASI  senses  the 
spinning  of  a  small  wind-driven  propeller  blade  using  a  cadmium  disulphide  optical  sensor. 
The  frequency  of  the  changes  in  light  intensity  sensed  is  transmitted  to  a  hand-held 
receiver  which  converts  the  frequency  to  airspeed  which  can  be  read  directly  in  real-time. 
The  manufacturer's  claimed  accuracy  is  +/-  0.5  feet  per  second.  The  airspeed  indicator 
was  installed  on  the  left  wing  of  the  Bluebird  by  mounting  it  on  the  end  of  a  boom  which 
extended  approximately  18  inches  (approximately  80  percent  of  c-bar)  in  front  of  the 
leading  edge  of  the  wing. 

Measuring  elevator  trim  was  done  in  an  indirect  manner.  The  flight  control 
transmitter  has  a  digital  display  which  can  show  elevator  trim  on  a  generic  scale  of -100  to 
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100.  Prior  to  flight  test,  elevator  deflection  in  degrees  was  measured  for  various  trim 
settings,  and  the  relationship  between  transmitter  displayed  trim  number  and  actual 
degrees  of  elevator  deflection  was  determined.  This  calibration  allowed  for  the 
determination  of  elevator  deflection  in  degrees  during  post  flight  analysis. 
C.      TEST  PROCEDURES 

In  order  to  fly  the  Bluebird  at  various  centers  of  gravity,  an  access  panel  in  the  aft 
fuselage  of  the  aircraft  was  modified  to  accept  added  weights.  The  aircraft  was  then 
configured  with  various  amounts  of  weight  and  the  location  of  the  center  of  gravity 
determined  for  each  configuration.  Centers  of  gravity  were  determined  by  weighing  each 
wheel  with  a  strain-gage  balance. 

A  total  of  seven  flights  were  flown  at  various  centers  of  gravity.  Flights  were  kept 
short  in  order  to  minimize  the  shift  in  the  center  of  gravity  due  to  fuel  burn.  The  location 
of  the  center  of  gravity  was  determined  both  before  and  after  the  flight,  and  the  average 
center  of  gravity  was  used  for  calculations.  Shifts  in  center  of  gravity  during  testing  were 
kept  to  one  percent  or  less  of  the  wing  chord. 

Each  flight  consisted  of  12  passes  at  various  airspeeds.  Airspeed  and  trim  setting 
for  each  pass  were  recorded  for  post  flight  analysis.  Additionally,  air  pressure,  air 
temperature,  and  aircraft  weight  were  noted  in  order  to  calculate  coefficient  of  lift. 

Post  flight  analysis  of  the  data  was  conducted  in  accordance  with  the  procedures 
described  in  Section  A,  Test  Description.  The  resulting  graphs  are  shown  in  Figures 
24-31. 
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Flight  1  of  7 
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Figure  25 
Flight  2  of  7 
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Figure  26 
Flight  3  of  7 
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Figure  27 
Flight  4  of  7 
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Figure  28 
Flight  5  of  7 
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Figure  29 
Flight  6  of  7 
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Flight  7  of  7 
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D.      RESULTS 

The  neutral  point  estimated  by  flight  test  was  found  to  be  about  54  percent  of  the 
mean  aerodynamic  chord.  The  neutral  point  estimated  through  conventional  calculations 
was  approximately  52.8  percent  of  the  mean  aerodynamic  chord. 

While  such  a  close  match  between  the  neutral  point  locations  found  by  each  method 
is  highly  encouraging,  one  must  not  overlook  the  scatter  in  the  data  presented  in  Figure 
3 1 .  Assuming  a  normal  distribution  of  the  data,  there  is  a  68  percent  confidence  that  the 
experimentally  determined  neutral  point  is  within  five  percent  of  the  wing  mean 
aerodynamic  chord  of  the  location  determined  experimentally.  Such  data  scatter  is 
unacceptable  for  accuracies  required,  and  the  tests  will  be  repeated  when  the  improved 
instrumentation  under  development  comes  on  line. 

For  the  flight  testing  done,  two  sources  of  potential  error  are  of  particular  interest. 
First,  it  was  very  difficult  to  ensure  a  perfectly  level  pass  of  the  aircraft  at  each  particular 
airspeed.  The  pilot  was  positioned  on  the  ground  and  the  plane  was  flying  approximately 
one  hundred  feet  above  him.  This  is  not  an  optimum  vantage  point  for  the  pilot.  Ideally, 
the  pilot  would  like  to  be  elevated  (such  as  in  a  tower)  such  that  the  aircraft  is  at  eye  level 
for  each  pass.  Future  plans  include  using  an  altitude-hold  autopilot  to  maintain  the  desired 
trim  conditions.  Also,  due  to  wind  gusts  and  possibly  small  unintentional  changes  in 
aircraft  attitude,  the  airspeed  was  not  always  steady,  and  therefore  somewhat  difficult  to 
read  consistently.  This  second  problem  will  be  overcome  by  the  improved  data  acquisition 
system. 
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VI.  SUMMARY 

A  determination  of  the  moments  of  inertia  by  the  component  contribution  method 
was  decided  to  be  too  cumbersome  for  analysis  of  the  complete  aircraft.  Moments  of 
inertia  were  found  through  a  compound  pendulum  analysis,  and  appear  reasonable.  Future 
changes  to  the  configuration  of  the  aircraft  can  be  accounted  for  by  considering  the 
contribution  of  the  component  added,  removed,  or  relocated. 

Initial  estimates  of  stability  and  control  derivatives  were  made  by  conventional 
aircraft-design-type  methods.  Such  methods  involve  considering  characteristics  of  the 
aircraft  such  as  lift-curve  slopes  of  airfoils  and  the  locations  of  particular  aircraft 
components  relative  to  each  other.  Programs  were  written  in  MATLAB  in  such  a  manner 
that  future  changes  to  the  aircraft  or  different  flight  conditions  can  be  accounted  for 
quickly  and  accurately.  The  initial  estimates  of  the  derivatives  can  be  used  in  the 
preliminary  design  of  a  flight  controller. 

The  initial  estimates  of  the  stability  and  control  derivatives  were  used  to  predict  the 
dynamic  response  of  the  aircraft  to  several  different  excitations.  MATLAB  programs 
were  written  to  predict  the  dynamics,  and  any  future  modifications  to  the  aircraft  or  flight 
conditions  can  be  accounted  for  easily.  Several  characteristics  of  the  aircraft's  responses 
are  worth  noting.  First,  the  short  period  longitudinal  response  was  found  to  be  heavily 
damped.  Little  or  no  evidence  of  the  short  period  response  can  be  observed  beyond 
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approximately  two  seconds.  Second,  the  long  period  longitudinal  response  was  found  to 
be  lightly  damped.  The  long  period  response  can  still  be  observed  relatively  easily  after 
100  seconds.  Similarly,  the  Dutch  roll  mode  is  also  fairly  lightly  damped  and  can  be 
observed  easily  for  ten  to  1 5  seconds.  Finally,  the  spiral  mode  was  found  to  be  divergent 
with  a  time  to  double  bank  angle  of  approximately  20  seconds.  While  this  mode  is  not 
particularly  dangerous  for  a  radio  controlled  aircraft  which  is  flown  visually  (open  loop), 
the  pilot  should  be  aware  of  this  tendency  of  the  aircraft  in  order  to  avoid  trim  problems. 
It  is  also  recommended  that  flight  controller  designs  incorporate  bank  angle  feedback  for 
wing  leveling. 

The  neutral  point  of  the  aircraft  was  found  to  be  53  to  54  percent  of  the  wing  mean 
aerodynamic  chord  by  two  different  methods.  The  data  collected  through  flight  test  were 
somewhat  scattered;  it  is  estimated  that  the  neutral  point  location  is  known  within  +/-  five 
percent  of  the  wing  mean  aerodynamic  chord  with  a  confidence  of  about  68  percent. 
Further  flight  testing  with  improved  instrumentation  is  recommended  to  raise  the 
confidence  of  the  neutral  point  estimation.  It  is  recommended  that  flight  tests  with 
improved  instrumentation  continue  to  verify  or  update  all  predicted  stability  and  control 
derivatives. 
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APPENDIX  A:  MATLAB  PROGRAMS  USED  TO  ESTIMATE 
STABILITY  AND  CONTROL  DERIVATIVES 


1.  blbrdfcl.m 


%  Eric  J.  Watkiss 

%  AA08 10  Thesis 

%  File  for  Bluebird  data  which  change  with  flight  condition 

%  blbrdfcl.m 

%  Last  Update:  02  MAR  94 


g  =  32.174; 
Wlmg  =  24.02 
Wrmg  =  23.91 
Wng=  11.07 

Umph  =  60 
Ufps  =  88.0 

rho  =  .002327 

Ixx=  12.58 
Iyy=  13.21 
Izz  =  19.99 
lxz  =  0 


%  Acceleration  due  to  gravity 
%  Weight  on  left  main  in  lbs 
%  Weight  on  right  main  in  lbs 
%  Weight  on  nose  gear  in  lbs 

%  Flight  speed  in  miles  per  hour 
%  Flight  speed  in  feet  per  second 

%  Air  density  in  slugs/(cubic  ft) 

%  Moment  of  inertia  about  x-axis 
%  Moment  of  inertia  about  y-axis 
%  Moment  of  inertia  about  z-axis 
%  Assumed!!!!!!!!!!!!!!!!!!!!!!! 


LD  =  8; 

thetanaut  =  0; 


%  Lift  to  drag  ratio 
%  Initial  pitch  angle 
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2.  bluebird. m 

%  Eric  J.  Watkiss 

%  AA08 10  Thesis 

%  File  for  Bluebird  data  which  are  fixed 

%  bluebird. m 

%  Last  Update:  02  MAR  94 


ac  =  .479; 

%  Aileron  chord  in  ft. 

ai  =  3; 

%  distance  from  centerline  to 

% 

inner  edge  of  aileron  in  ft. 

alpha01  =  -6.5*pi/180; 

%  a.o.a.  for  zero  lift  (radians)ao  =  6, 

% 

distance  from  centerline  to 

% 

outer  edge  of  aileron  in  ft. 

b=  12.42; 

%  Span  of  wing  in  ft 

bt  =  3.67; 

%  Span  of  horizontal  tail  in  ft. 

bv=  1.208; 

%  Height  of  vertical  tail  in  ft. 

cbar=  1.802; 

%  Mean  aerodynamic  chord  (mac.) 

% 

in  feet 

CLalphaaf\v  =  5.443; 

%  Lift  curve  slope  of  wing 

% 

airfoil  (GO  769)  in  per 

% 

radian 

CLalphaaft  =  5.587; 

%  Lift  curve  slope  of  horizontal 

% 

tail  airfoil  (NACA  4412)  in  per 

% 

radian 

CLalphaafV  =  2*pi; 

%  Lift  curve  slope  of  vertical 

% 

tail  airfoil  (flat  plate)  in  per 

% 

radian 

CMac  =  -.06; 

%  Coefficient  of  moment  about 

% 

aero.  ctr.  (GO  769) 

ct  =  1  281; 

%  mac.  of  horizontal  tail  in  ft. 

c4tail  =  7.878; 

%  Location  of  quarter  chord  of 

% 

horizontal  tail  in  feet  from 

% 

firewall 

c4wing  =  2.497; 

%  Location  of  quarter  chord  of 

% 

wing  in  feet  from  firewall 

daOdde  =  .625; 

%  Section  flap  effectiveness 

% 

for  33%  flap  (elevator) 

% 

Abbott  and  Doenhoffp.  190 

daOddr  =  .675; 

%  Section  flap  effectiveness 

% 

for  38%  flap  (rudder) 

% 

Abbott  and  Doenhoffp.  190 

deda  =  .4; 

%  Downwash  angle  derivative 
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% 

estimated  from  Perkins/Hage 

Df=  1; 

%  Depth  of  fuse,  in  ft. 

eO  =  0; 

%  Assumed  epsilon  naught 

ee=8; 

%  Assumed  span  efficiency  factor 

g  =  32.174; 

%  gravitational  constant 

hac  =  .245; 

%  Location  in  percent  chord  of 

% 

aero.  ctr.  (NACA4412) 

it  =  4.83*pi/180; 

%  Incidence  angle  of  hor.  tail 

lewing  =  2.047; 

%  Location  of  leading  edge  of  wing 

% 

in  feet  from  firewall 

letail  =  7.557, 

%  Location  of  leading  edge  of 

% 

horizontal  tail  in  feet  from 

% 

firewall 

mg  =  37.595/12; 

%  Location  of  main  gear  in  ft 

% 

from  firewall 

ng=.  75/12; 

%  Location  of  nose  gear  in  ft 

% 

from  firewall 

S  =  22.380; 

%  Reference  (wing)  area  in  sq.  ft. 

Sr=.547; 

%  Rudder  area  in  sq.  ft. 

St  =  4.701; 

%  Horizontal  tail  area  in  sq.  ft. 

Sv  =  1.277; 

%  Vertical  tail  area  in  sq.  ft. 

Wf  =  .67; 

%  Width  of  fuse,  in  ft. 

ybar  =  b/4; 

%  Spanwise  location  of  max. 

zv  =  .5; 

%  Vert,  tail  height  to  m.a.c. 

% 

(estimated) 

Zwf=.5; 

%  Verticle  height  of  wing 

% 

above  fuse.  C.L.  in  ft. 
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3.  dderiv.m 

%  Eric  J.  Watkiss 

%  AA08 10  Thesis 

%  File  to  calculate  dimensional  derivatives 

%  dderiv.m 

%  Last  Update:   12  FEB  94 

%  Run  nondimensional  derivative  program 
ndderiv 

%  Calculate  dynamic  pressure 

qbar  =  5*rho*UfpsA2;  %  ft  lbs 

Malpha  =  CMalpha*qbar*S*cbar/Iyy;  %  per  secondA2 

Mq  =  CMq*(cbar/(2*Ufps))*qbar*S*cbar/Iyy; 

%  per  second 

Malphadot  =  CMalphadot*(cbar/(2*Ufps))*qbar*S*cbar/Iyy, 

%  per  second 

Xu  =  -2*CD*qbar*S/(m*Ufps);  %  per  second 

Zu  =  -2*CL*qbar*S/(m*Ufps);  %  per  second 

Zalphadot  =  CLalphadot*(cbar/(2*Ufps))*(qbar*S/m); 

%  ft  per  second 

Zq  =  CLq*(cbar/(2*Ufps))*(qbar*S/m);  %  ft  per  second 

Mu  =  0;  %  per  ft  second 

Xde  =  - 1  *CDde*qbar*S/m;  %  ft  per  secondA2 

Zde  =  - 1  *CLdelta*qbar* S/m;  %  ft  per  secondA2 

Mde  =  CMde*qbar*S*cbar/Iyy;  %  per  secondA2 

Xalpha  =  (CL  -  CDalpha)*qbar*  S/m;  %  ft  per  secondA2 

Zalpha  =  - 1  *(CLalphaw+CD)*qbar*S/m;  %  ft  per  secondA2 
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YB 

=  CyB*qbar*S/m; 

%  ft/secondA2 

LB  = 

=  ClB*qbar*S*b/Ixx; 

%  l/secondA2 

NB 

=  CnB*qbar*S*bAzz; 

%  l/secondA2 

Yp  = 

=  Cyp*b*qbar*S/(2*Ufps 

*m); 

%  ft/sec 

Yr  = 

Cyr*b*qbar*S/(2*UfpsI( 

:m); 

%  ft/sec 

Lp  = 

=  Clp*(b/(2*Ufps))*qbar* 

S*b/Ixx; 

%  1/sec 

Np  = 

=  Cnp*(b/(2*Ufps))*qbar 

*S*b/Izz; 

%  1/sec 

Lr  = 

Clr*(b/(2*Ufps))*qbar*< 

S*b/Ixx; 

%  1/sec 

Nr  = 

:  Cnr*(b/(2*Ufps))*qbar* 

S*b/Izz; 

%  1/sec 

Ydr 

=  -l*Cydr*qbar*S/m; 

%  ft/secA2 

Yda 

=  0; 

%  ft/secA2 

Ldr: 

=  Cldr*qbar*S*b/Ixx; 

%  l/secA2 

Lda 

=  Clda*qbar*S*b/Ixx; 

%  l/secA2 

Ndr 

=  Cndr*qbar*S*b/Izz; 

%  l/secA2 

Nda 

=  Cnda*qbar*S*b/Izz; 

%  l/secA2 

Malphaprime  =  Malpha  +  Malphadot*(Zalpha/Ufps); 

Mqprime  =  Mq  +  Malphadot; 

Mdeprime  =  Mde  +  Malphadot  *(Zde/Ufps); 
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4.  godallas.m 

%  Eric  J.  Watkiss 

%  AA08 10  Thesis 

%  File  to  get  values  of  dimensional  and  nondimensional  derivatives 

%  godallas.m 

%  Last  Update:  04  FEB  94 

!del  bluebird. dia 
diary  bluebird. dia 
diary  on 

%  Run  programs  to  calculate  derivatives 
dderiv 

%  Nondimensional  Derivatives 

CL 

CD 

CDO 

CLalphaw 

CMalpha 

CDalpha 

CLalphadot 

CMalphadot 

CLq 

CMq 

CLdelta 

CDde 

CMde 

CyB 

CnB 

C1B 

Cyp 

Cnp 

Clp 

Cyr 

Cnr 

Clr 

Cyda 

Cnda 

Clda 

Cydr 

Cndr 


49 


Cldr 
CDq 

Cyq 

%  Longitudinal  Dimensional  Derivatives 

Xu 

Xalpha 

Xde 

Zu 

Zalpha 

Zalphadot 

Zq 

Zde 

Mu 

Malpha 

Malphadot 

Mq 

Mde 

%  Lateral/Directional  Dimensional  Derivatives 
YB 

Yp 

Yr 

Yda 

Ydr 

LB 

Lp 

Lr 

Lda 

Ldr 

NB 

Np 

Nr 

Nda 

Ndr 

diary  off 
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5.  ndderiv.m 

%  Eric  J.  Watkiss 

%  AA08 10  Thesis 

%  File  to  calculate  nondimensional  derivatives 

%  ndderiv.m 

%  Last  Update:   12  FEB  94 

%  Load  Bluebird  data  with  flight  condition 
physcalc 

%  Calculate  coefficients  of  lift  and  drag 
CL  =  W/(.5*rho*UfpsA2*S); 
CD  =  CL/LD; 

%  Calculate  lift  curve  slope  of  wing  in  per  radian 
CLalphaw  =  CLalphaafw/(l+CLalphaafw/(pi*ee*AR)); 

%  Calculate  lift  curve  slope  of  horizontal  tail  in  per  radian 
CLalphat  =  CLalphaaft/(l+CLalphaaft/(pi*ee*ARt)); 

%  Calculate  lift  curve  slope  of  vertical  tail  in  per  radian 
CLalphav  =  CLalphaafV/(l+CLalphaafv/(pi*ee*ARv)); 

%  Calculate  change  in  nor.  tail  lift  with  change  in  elevator 
dcLtdde  =  daOdde  *  CLalphat;  %  per  radian 

%  Calculate  change  in  vert,  tail  lift  with  change  in  rudder 
dcLvddr  =  daOddr  *  CLalphav;  %  per  radian 

%  Calculate  zero  lift  pitching  moment 
CMO  =  CMac  +  VH  *  CLalphat  *  (it  +  eO); 

%  Calculate  CMalpha  in  per  radian 

CMalpha  =  CLalphaw*((h-hac)-VH*(CLalphat/CLalphaw)*(l-deda)); 

%  Calculate  change  in  aircraft  lift  with  change  in  elevator 
CLdelta  =  dcLtdde*(St/S);      %  per  radian 

%  Calculate  chng  in  aircraft  pitching  moment  w.  chng  in  elevator 
CMde  =  - 1  * VH*dcLtdde;       %  per  radian 

%  Calculate  angle  of  attack  and  elevator  angle  for  trimmed  flight 
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% 

%         CM  =  CMO  +  CMalph 

ia*alpha  +  CMde*de 

%         CI  =  CLalphaw*alpha 

+  CLdelta*de 

% 

% 

%         |  CLalphaw    CLdelta 

|  |  alpha  |      |  CL        | 

%         | 

II            1  =  1 

%         |  CMalpha      CMde 

||      de    |      |  -CMO  | 

% 

%                   A            *      X 

=     C 

% 

A  =  [  CLalphaw  CLdelta 

CMalpha     CMde    ]; 

C  =  [     CL 

-l*CMO]; 

X  =  inv(A)*C; 

atrim  =  X(l,l); 

%  trim  a.o.a.  in  radians 

etrim  =  X(2,l); 

%  trim  elevator  in  radians 

%  Calculate  change  in  yawing  moment  with  change  in  rudder 

%  "rudder  power" 

%  assumes  VF/Vinfinity  =  1 

Cndr  =  -1  *  W*dcLvddr;        %  in  per  radian 

%  Calculate  CnB  contribution  from  vert,  tail 

%  CnB  =  CLalphav*W*(VF/Vinfinity)A2*(l-dsigma/dbeta) 

%  assumes  VF/Vinfinity  =  1  and  dsigma/dbeta  =  0 

CnB  =  CLalphav*  W;%  in  per  radian 

%  Calculate  change  in  rolling  moment  with  change  in  sideslip 

%  First  calculate  dihedral  contribution  from  wing 
%  Raymer  p.  439 

ClBwf  =  -1.2*sqrt(AR)*Zwf*(Df+Wf)/bA2; 
ClBw  =  ClBwCL*CL+ClBwf; 

%  Next  calculate  contribution  from  fin 

%  ClBv  =  -l*Clalphav*Vprime*(VF/Vinfinity)A2*(l -dsigma/dbeta) 

%  Assume  VF/Vinfinity  =  1  and  dsigma/dbeta  =  0 

ClBv  =  - 1  *CLalphav* Vprime;  %  in  per  rad 


%  Combine  ClBg  and  ClBv  into  C1B 
C1B  =  ClBw  +  ClBv;   %  in  per  rad 
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%  Calculate  "aileron  power",  Clda 
%  See  Smetana  pp.  139-141 
Cldatau  =  Cldatauo  -  Cldataui; 
Clda  =  Cldatau*tau;    %  in  per  radian 

%  Calculate  change  in  yawing  moment  w.  aileron  deflection 
Cnda  =  2*K*CL*Clda;  %  in  per  radian 

%  Calculate  side  force  due  to  yaw 

%  By  Smetana  p.  107 

CyB  =  -.31,      %  in  per  radian 

%  Calculate  side  force  due  to  rudder 

Cydr  =  CLalphav*taur*Sv/S;  %  in  per  radian 

%  Calculate  side  force  due  to  aileron 
%  By  Smetana,  p.  138 
Cyda  =  0; 

%  Calculate  rolling  moment  due  to  rudder 
Cldr  =  Cydr*zv/b;        %  in  per  radian 

%  Calculate  change  in  drag  due  to  change  in  elevator 

%  Smetana  pp.  95-100 

%  Using  Figure  26  at  8  degrees  aoa 

CDde  =  ((.155-.047)/(20*pi/180))*St/S;        %  in  per  radian 

%  Calculate  change  in  drag  with  change  in  aoa 

%  Smetana  pp.  64-65 

%  Assuming  dCDO/dalpha  is  negligible 

CDalpha  =  2*CL*CLalphaw/(pi*ee*AR);      %  in  per  radian 

%  Calculate  change  in  pitching  moment  w.r.t.  alphadot 
%  Smetana  pp.  78-81,  etat  assumed  =  1 
CMalphadot  =  -2*CLalphat*deda*(ltprime/cbar)*... 
(lt/cbar)*(St/S);  %  in  per  radian 

%  Calculate  change  in  lift  with  pitch  rate 

%  Smetana  pp.  82-85 

%  Neglecting  wing  contribution,  assuming  etat  =  1 

CLq  =  2*(lt/cbar)*CLalphat*(St/S);  %  in  per  radian 
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%  Calculate  change  in  lift  with  alphadot 

%  Smetana  pp.  75-76 

CLalphadot  =  -1  *CMalphadot/(lt/cbar);         %  in  per  radian 

%  Calculate  change  in  pitching  moment  w.  pitch  rate 

%  Smetana  pp.  87-88 

%  Assuming  etat  =  1 

CMq  =  -2*(cbar/4-h)*abs(cbar/4-h)*CLalphaw/(cbarA2)  - ... 

2*(lt/cbar)A2*CLalphat*(St/S);  %  in  per  radian 

%  Calculate  roll  damping 

%  Smetana  pp.  122-125 

%  Neglecting  contribution  from  vertical  tail 

Clp  =  -.475*(AR+4)/(2*pi*AR/CLalphaw+4);  %  in  per  radian 

%  Calculate  change  in  yawing  moment  due  to  rolling 

%  Smetana  pp.  126-129 

%  Neglecting  contribution  from  vertical  tail 

Cnp  =  -1  *CL/8,  %  in  per  radian 

%  Calculate  change  in  side  force  with  yaw  rate 

%  From  Schmidt  p.  3-23 

%  Assume  etat  =  1 

Cyr  =  2*  W*CLalphav;  %  in  per  radian 

%  Calculate  change  in  rolling  moment  w.  yaw  rate 

%  Schmidt  p.  3-24 

%  Tail  contribution 

Clrv  =  (zv/b)*Cyr;       %  in  per  radian 

%  Wing  contribution 

Clrw  =  CL/4;   %  in  per  radian 

%  Total 

Clr  =  Clrv  +  Clrw;       %  in  per  radian 

%  Calculate  yaw  damping 

%  Schmidt  p.  3-25 

%  Tail  contribution 

Cnrv  =  -l*(lv/b)*Cyr;  %  in  per  radian 

%  Wing  contribution  from  Smetana  p.  136 

CDO  =  CD-CLA2/(pi*ee*AR); 

Cnrw  =  -02*CLA2-.3*CD0;  %  in  per  radian 

%  Total 

Cnr  =  Cnrv  +  Cnrw;    %  in  per  radian 
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%  The  following  3  derivatives  are  negligible  and  taken  to  be  0 
CDq  =  0;         %  in  per  radian 
Cyq  =  0;  %  in  per  radian 

Cyp  =  0;  %  in  per  radian 

%  A  few  misc.  calculations 

%  Static  Margin/Neutral  Point 
statmar  =  CMalpha/(-l*CLalphaw) 
hn  =  statmar  +  h 
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6.   physcalc.m 

%  Eric  Watkiss 

%AA08 10  Thesis 

%  File  to  calculate  physical  considerations 

%  physcalc.m 

%  Last  Update:  04  FEB  94 

%  Load  fixed  Bluebird  data 
bluebird 

%  Load  flight  condition 
blbrdfcl 

%  Calculate  aircraft  weight 
W  =  Wlmg  +  Wrmg  +  Wng, 

%  Calculate  aircraft  mass 
m  =  W/g; 

%  Calculate  aspect  ratio  of  wing 
AR  =  bA2/S; 

%  Calculate  aspect  ratio  of  hor.  tail 
ARt  =  btA2/St; 

%  Calculate  aspect  ratio  of  vert,  tail 
ARv  =  bvA2/Sv; 

%  Calculate  longitudinal  center  of  gravity 

h  =  ((ng*Wng  +  mg*(Wlmg+Wrmg))AV-lewing)/cbar; 

%  Calculate  "tail  length"  from  e.g.  to  horizontal  tail  a.c. 
%  same  for  horizontal  and  vertical 
It  =  c4tail  -  (lewing  +  h*cbar); 
lv  =  It; 

%  Calculate  "tail  length"  from  c/4  wing  to  c/4  tail 
ltprime  =  c4tail  -  c4wing; 

%  Calculate  hor.  tail  volume  coefficient 
VH  =  lt*St/(S*cbar); 
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%  Calculate  vert,  tail  volume  coefficient  (yaw) 
W  =  lv*Sv/(b*S); 

%  Calculate  vert,  tail  volume  coefficient  (roll) 
Vprime  =  zv*Sv/(b*S); 

%  Unit  antisymmetrical  angle  of  attack  for  outer  and  inner 

%  edge  of  aileron  (See  Smetana  p.  141) 

antisymo  =  ao/(b/2); 

Cldatauo  =  .74; 

antisymi  =  ai/(b/2); 

Cldataui  =  .23; 

cacw  =  ac/cbar; 

tau  =  .52; 

%  for  yawing  moment  due  to  aileron,  see  p.  142,  Smetana 
eta  =  ai/(b/2); 
K  =  -.17; 

%  for  side  force  due  to  rudder  deflection,  see  Smetana  p.  145 
vratio  =  Sr/Sv; 
taur  =  .62; 

%  for  rolling  moment  due  to  sideslip,  See  Raymer,  Fig.  16.21,  p.  439 
ClBwCL  =  -.04; 
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APPENDIX  B:  MOMENT  OF  INERTIA  CALCULATIONS 

For  the  formula 

r    _     Wm  +  2<o(L  short  +  Llong)¥m+s  n  wmZm  2(Dfr2  r  2  "^ 

1    -  4„2  'M+S  g  3g  \L  SHORT +L  LONG  J 

the  following  data  were  used  for  extraction  of  the  moments  of  inertia. 

Fixed  values: 

Weight  of  the  model,  WM  =  58.45  lbs 

Weight  per  unit  length  of  chain,  co  =  0.06133  lbs/ft 

Length  of  short  chain,  Lshort  =  13.00ft 

Length  of  long  chain,  Llong  =  15.00ft 

Gravitational  constant,  g  =  32.1472  ft/sec2 
(adjusted  for  latitude  and  elevation) 

Variable  values: 

L^:      Distance  from  pivot  to  center  of  gravity 

of  model  and  support,  Zm+s  =  12.95ft 

Distance  from  pivot  to  center  of  gravity 

of  model,  Zm  -  13.33  ft 

Average  period  of  model  and  support,  Pm+s  -  4.109  sec 

L,:      Distance  from  pivot  to  center  of  gravity 

of  model  and  support,  ZM+s  -  12.01  ft 

Distance  from  pivot  to  center  of  gravity 

of  model,  Zm  =  12.35  ft 

Average  period  of  model  and  support,  Pm+s  =  3.976  sec 

1^:      Distance  from  pivot  to  center  of  gravity 

of  model  and  support,  ZM+S  =  12.49ft 

Distance  from  pivot  to  center  of  gravity 

ofmodel,  Zm  -  12.81  ft 

Average  period  ofmodel  and  support,  Pm+s  =  4.077  sec 
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APPENDIX  C:  MATLAB  PROGRAMS  USED  TO  ESTIMATE 
AIRCRAFT  DYNAMICS 


A.      LONGITUDINAL  DYNAMICS 

1.       longnat.m 

%  Determines  the  short  and  long  period  natural  response 
%  using  the  full  longitudinal  plant 

clear 

dderiv 

!del  longnat.dia 

diary  longnat.dia 

%  Inertial  Matrix 

in  =  [Ufps  0  0    0; 

0        Ufps-Zalphadot  0    0; 

0        -l*Malphadot    1    0; 

0  0  0    1]; 

%  Aircraft  Matrix 

an  =  [Ufps*Xu     Xalpha       0  -l*g*cos(thetanaut); 

Ufps*Zu    Zalpha    Ufps+Zq  -l*g*sin(thetanaut); 

Ufps*Mu     Malpha     Mq  0; 

0                     0             1  0]; 

a  =  inv(in)*an; 

b  =[1;1;1;1]; 

c  =  eye(size(a)); 

d  =[0;0;0;0]; 

t  =  0:01:100; 

[y,x,t]  =impulse(a,b,c,d,l,t); 

clg 

figure(l) 

%plotting  alpha  and  q,  short  period 
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plot(t,y(:,2));axis([0  15  -1.5  1.5]);grid;gtext('alpha') 

hold  on 

plot(t,y(:,3));gtext('q') 

xlabel(Time,  sec');ylabel('Response  Gain') 

VotitleCNatural  Response') 

hold  off 

pause 

figure(2) 

%plotting  u/U  and  theta,  phugoid 

plot(t,y(:,l));axis([0  100  -2  2]);grid;gtext('u/U') 

hold  on 

plot(t,y(:,4));gtext('theta') 

xlabel('Time,  sec'),ylabel('Response  Gain') 

%title(rNatural  Response'); 

hold  off 

damp(eig(a)) 

diary  off 
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2.       stepper,  m 

%  Determines  the  short  and  long  period  step  response 
%  using  the  full  longitudinal  plant 

clear 
dderiv 

%  Inertial  Matrix 

in  =  [Ufps         0  0    0; 

0    Ufps-Zalphadot  0    0; 

0     -l*Malphadot    1    0; 

0  0  0    1]; 

%  Aircraft  Matrix 

an  =  [Ufps*Xu     Xalpha    0       -l*g*cos(thetanaut); 

Ufps*Zu     Zalpha    Ufps+Zq  -l*g*sin(thetanaut); 

Ufps*Mu     Malpha    Mq  0; 

0  0        1  0]; 

%  Control  Matrix 

bn  =  [Xde  Zde  Mde  0]'; 

a  =  inv(in)*an; 

b  =  inv(in)*bn; 

c  =  eye(size(a)); 

d  =[0  0  0  0]'; 

t  =  0:01:100; 

[y,x,t]  =step(a,b,c,d,l,t); 

clg 

figure(l) 

%plotting  alpha  and  q,  short  period 

plot(t,y(:,2)) 

axis([0  15-5  3]) 

gridigtextCalpha') 

hold  on 

plot^yCS^gtextCq') 

xlabel('Time,  sec^iylabelCResponse  Gain') 

%title('Short  Period  Step  Response') 

hold  off 
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pause 

figure(2) 

%plotting  u/U  and  theta,  phugoid 

plot(t,y(:,l)) 

axis 

grid;gtext('u/U') 

hold  on 

plot(t,y( :  ,4));gtext('theta') 

xlabel(Time,  sec'),ylabel('Response  Gain') 

%title('Long  Period  Step  Response'); 

hold  off 

pause 

figure(3) 

%plotting  alpha  and  q,  long  period 

plot(t,y(:,2)) 

grid;gtext('alpha') 

hold  on 

plot(t,y(:,3));gtext('q') 

xlabel(Time,  sec');ylabel('Response  Gain') 

%title('Long  Period  Step  Response') 

hold  off 
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3.       nstep.m 

%    This  program  uses  STEP  function  to  determine 
%    normal  acceleration  response  to  a  unit  step  input 

clear 
dderiv 

%  Inertial  Matrix 

in  =  [Ufps         0  0    0; 

0    Ufps-Zalphadot  0    0; 

0     -l*Malphadot    1    0; 

0  0  0    1]; 

%  Aircraft  Matrix 

an  =  [Ufps*Xu     Xalpha    0        -l*g*cos(thetanaut); 

Ufps*Zu     Zalpha    Ufps+Zq  -l*g*sin(thetanaut); 

Ufps*Mu     Malpha     Mq  0; 

0  0        1  0]; 

%  Control  Matrix 

bn  =  [Xde  Zde  Mde  0]'; 

a  =  inv(in)*an; 

b  =  inv(in)*bn; 

c  =  [  Zu  Zalpha  Zq  0  ]; 

d  =Zde; 

t    =0:0.05:15; 

[y,x,t]  =step(a,b,c,d,l,t); 

fori=  l:length(t) 

y2  =y/32.174*pi/180; 
end 

%  plotting  g's/deg 

plot(t,y2);grid 

xlabel(Time,  sec');ylabelCNormal  Acceleration:  g  s/deg') 

%title(' Acceleration  Response  to  Unit  Step') 
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4.       homogen.m 

%  Determines  the  short  and  long  period  homogeneous  response 
%  using  the  full  longitudinal  plant 

clear 
dderiv 

%  Inertial  Matrix 

in  =  [Ufps         0  0    0; 

0    Ufps-Zalphadot  0    0; 

0     -l*Malphadot    1    0, 

0  0  0    1]; 

%  Aircraft  Matrix 

an  =  [Ufps*Xu     Xalpha     0        -l*g*cos(thetanaut); 

Ufps*Zu     Zalpha    Ufps+Zq  -l*g*sin(thetanaut); 

Ufps*Mu     Malpha     Mq  0, 

0  0        1  0]; 

%  Control  Matrix 

bn  =  [Xde  Zde  Mde  0]'; 

a  =  inv(in)*an; 

b  =  inv(in)*bn; 

c  =  eye(size(a)); 

d  =[0  0  0  0]'; 

t  =  0:01:100; 

[y,x,t]  =step(a,b,c,d,l,t); 

xO  =y(find(t=15),:); 

xO  =  xO/xO(2)*5/57.3 
[y,x,t]  =  initial(a,b,c,d,x0,t); 

clg 

figure(l) 

%plotting  alpha  and  theta,  short  period 

plot(t,y(:,2)) 

%title('Short  Period  Homogeneous  Response'); 

grid 

axis([0  15-.2.2]); 
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gtext('alpha  -  rad'); 

hold  on 

plot(t,y(:,3));gtext('q  -  rad/sec'); 

xlabel(Time,  sec1);  ylabel('Response  Gain'); 

hold  off 

pause 

figure(2) 

%plotting  u/U  and  theta,  long  period 

plot(t,y(:,l)) 

%title('Long  Period  Homogeneous  Response'); 

grid 

axis; 

gtextCu/U'); 

hold  on 

plot(t,y(:,4));gtext('theta  -  rad'); 

xlabel(Time,  sec');  ylabel('Response  Gain'); 

hold  off 
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5.       doublet,  m 

%    This  program  uses  EXPM  function  to  find  response  to  pitch  doublet 

clear 
dderiv 

%  Inertial  Matrix 

in  =  [Ufps         0  0    0; 

0    Ufps-Zalphadot  0    0, 

0     -l*Malphadot    1    0; 

0  0  0    1]; 


%  Aircraft  Matrix 

an  =  [Ufps*Xu     Xalpha     0       -l*g*cos(thetanaut); 

Ufps*Zu     Zalpha    Ufps+Zq  -l*g*sin(thetanaut); 

Ufps*Mu     Malpha     Mq  0; 

0  0        1  0]; 

bn  =  [  Xde  Zde  Mde  0  ]'; 

a    =  inv(in)*an, 
b    =  inv(in)*bn, 

tl  =1.0;     %  start  of  doublet,  sec 
t2  =  2.0;     %  midpoint  of  doublet,  sec 
t3  =3.0;     %  end  of  doublet,  sec 

dO  =  -1;  %unitelevat<     nput  (1  rad)  [t.e.u] 

tim  =15;  %  set  end  time 

t    =0:0.05:tim; 

x  =  zeros(4,length(t));      %  initialize  solution  matrices 
xl  =  zeros(4,length(t)); 
x2  =  zeros(4,length(t)); 
x3  =  zeros(4,length(t)); 

fori=  l:length(t) 

ift(i)>=tl 
de(i)    =d0; 
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xl(:,i)  =  dO*(expm(a*(t(i)-tl))  -  eye(size(a)))*inv(a)*b; 
end 

ift(i)>=t2 

de(i)    =  de(i)-2*d0; 

x2(:,i)  =  -2*d0*(expm(a*(t(i)-t2))  -  eye(size(a)))*inv(a)*b; 
end 

ift(i)>=t3 

de(i)    =de(i)  +  dO; 

x3(:,i)  =  dO*(expm(a*(t(i)-t3))  -  eye(size(a)))*inv(a)*b; 
end 

x(:,i)  =xl(:,i)  +  x2(:,i)  +  x3(:,i); 

end 

ymin  =  -2*abs(d0); 
ymax  =  2*abs(d0); 
V  =  [0  tim    ymin    ymax]1; 

figure(l) 

plot(t,de);grid;axis(V) 

xlabel(Time,  sec');  ylabel('Elevator,  Rad'); 

%title('Elevator  Input') 

pause;  axis 

figure(2) 

%plotting  alpha  and  pitch  rate 

plot(t,x(2,:)) 

%title('Doublet  Response') 

grid;gtext('alpha') 

hold  on 

plot(t,x(3,:));gtext(,q') 

xlabel('Time,  sec');ylabelCResponse  Gain') 

hold  off 
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sp_bode.m 


%    This  program  uses  the  BODE  function  to  find  the 

%    short-period  response  (transfer-function  gain  and  phase) 

%    to  a  harmonic  input 

clear 
dderiv 

%  Inertial  Matrix 

in  =  [Ufps         0  0    0; 

0    Ufps-Zalphadot  0    0; 

0     -l*Malphadot    1    0; 

0  0  0    1]; 


%  Aircraft  Matrix 

an  =  [Ufps*Xu     Xalpha     0       -l*g*cos(thetanaut); 

Ufps*Zu     Zalpha    Ufps+Zq  -l*g*sin(thetanaut); 

Ufps*Mu     Malpha     Mq  0; 

0  0        1  0]; 

bn  =  [  Xde  Zde  Mde  0  ]'; 

a  =  inv(in)*an; 

b  =  inv(in)*bn; 

c  =  eye(size(a)); 

d  =  zeros(size(b)); 

w  =  logspace(-l,2,150); 
[mag,phase,w]  =  bode(a,b,c,d,l,w); 

for  i  =  l:length(w) 
forj  =  l:4 
ifphase(i,j)  >0 

phase(ij)  =  phase(i,j)  -  360; 
end 
end 
end 

clg 

figure(l) 

%  Plotting  alpha  and  pitch-rate  gains 
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semilogx(w,mag(:,l));grid;gtext(,alpha,) 

%title('Short-Period  Frequency  Response') 

hold  on 

semilogx(w,mag( :  ,2));gtext('q') 

xlabel('Frequency,  rad/sec');ylabel('Gairi) 

hold  off 

pause 

figure(2) 

%  Plotting  alpha  and  pitch-rate  phase  angles 

semilogx(w,phase(:,l));grid;gtext('alpha') 

%title('Short-Period  Frequency  Response') 

hold  on 

semilogx(w,phase(:,2));gtext('q') 

xlabel('Frequency,  rad/sec');ylabel('Phase,  deg') 

hold  off 
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7.       n_bode.m 

%    This  program  uses  the  BODE  function  to  find 

%    the  normal-acceleration  response 

%    (transfer-function  gain  and  phase)  to  a  harmonic  input 

clear 
dderiv 

%  Inertial  Matrix 

in  =  [Ufps         0  0    0; 

0    Ufps-Zalphadot  0    0; 

0     -l*Malphadot    1    0; 

0  0  0    1]; 

%  Aircraft  Matrix 

an  =  [Ufps*Xu     Xalpha     0        -l*g*cos(thetanaut); 

Ufps*Zu     Zalpha    Ufps+Zq  -l*g*sin(thetanaut); 

Ufps*Mu     Malpha     Mq  0; 

0  0        1  0]; 

%  Control  Matrix 

bn  =  [Xde  Zde  Mde  0]'; 

a  =  inv(in)*an, 

b  =  inv(in)*bn; 

c  =  [  Zu  Zalpha  Zq  0  ]; 

d  =Zde; 

w  =  logspace(-l,2,150); 
[mag,phase,w]  =  bode(a,b,c,d,l,w); 

for  i  =  l:length(w) 

mag(i)    =  mag(i)/g*pi/180;  %  converting  to  g's/deg 

if  phase(i)  >  0 
phase(i)  =  phase(i)  -  360; 

end 
end 

%  Plotting  load  factor/deg  elevator  input 

figure(l) 

semilogx(w,mag);grid 

xlabel(Trequency,  rad/sec');  ylabelCNormal  Acceleration  Gain') 
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%titleCNormal  Acceleration  Response  to  Harmonic  Input') 
pause 

%  Plotting  normal  acceleration  phase  angle 

figure(2) 

semilogx(w,phase);grid 

xlabel('Frequency,  rad/sec');  ylabelCNormal  Acceleration  Phase,  deg') 

%title(TSformal  Acceleration  Response  to  Harmonic  Input') 
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B.       LATERAL-DIRECTIONAL  DYNAMICS 

1.       lat_dir.m 

%  Solves  the  full  4x4  lateral-directional  response 
clear 
dderiv 

!del  latdir.dia 
diary  latdir.dia 
in  =  [Ufps      0         0        0 
0       1         0    -l*Ixz/Ixx 
0       0         10 
0  -l*Ixz/Izz    0        1]; 
an  =  [YB     Yp    g*cos(thetanaut)     Yr-Ufps 

LB     Lp  0  Lr 

0       10  0 

NB     Np  0  Nr]; 

a  =  inv(in)*an; 
[wn,zeta]  =  damp(a) 
[x,d]  =  eig(a) 
r  =  [0  0  0  0]; 

for  i  =  1 :4 

r(i)  =  d(i,i); 
end 

x2  =  zeros(4,2); 

forj  =  1:4 

x2(j,l)  =  abs(x(j,l)); 

x2G,2)=180/pi*angle(xO,l)); 
end 

xx  =  x2(l,l); 
xy  =  x(3,2); 
xz  =  x(4,3); 


x3( 
x3( 
x3( 
x3( 
x3 


,l)  =  x2(:,l)/xx; 
,2)  =  x2(:,2)-x2(l,2); 
,3)  =  x(:,3)/xy; 
,4)  =  x(:,4)/xz; 


diary  off 
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2.       rudkick.m 

%  Response  to  unit  rudder  impulse 

clear 
dderiv 

in  =  [Ufps      0         0        0 
0       1         0    -l*Ixz/Ixx 
0       0         10 
0  -l*Ixz/Izz    0        1]; 

%  Plant  Matrix 

an  =  [YB     Yp    g*cos(thetanaut)     Yr-Ufps 
LB     Lp  0  Lr 

0      10  0 

NB     Np  0  Nr]; 

%  Control  Matrix 
bn  =  [Ydr  Yda 

Ldr  Lda 
0     0 

Ndr  Nda]; 

a  =  inv(in)*an; 
b  =  inv(in)*an; 
c  =  eye(size(a)); 
d  =  zeros(size(b)); 

t  =0:.05:10; 

[x,y,t]  =  impulse(a,b,c,d,l,t); 

%  plotting  sideslip  (beta)  and  bank  angle  (phi) 

plot(t,y(:,  l));grid;gtext('beta') 

hold  on 

plot^yCS^gtextCphi') 

xlabel(Time,  sec'XylabelCResponse  Gain1) 

%title('Rudder  Kick  Response') 

hold  off 

pause 

%  integrate  p  for  final  bank  angle 
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phi  =  0; 

t    =0:.05:100; 

[x,y,t]  =  impulse(a,b,c,d,l,t); 

fori=  llength(t) 

phi  =  phi  +  y(i,2)*05, 
end 

phi 
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3.       sideslip,  m 

%  Solves  the  full  4x4  lateral-directional  response  for 
%  initial  sideslip  condition 

clear 
dderiv 

in  =  [Ufps      0         0        0 
0       1         0    -l*Ixz/Ixx 
0       0         10 
0  -l*Ixz/Izz    0        1]; 

an  =  [YB     Yp    g*cos(thetanaut)     Yr-Ufps 
LB     Lp  0  Lr 

0      10  0 

NB     Np  0  Nr]; 

bn  =  [Ydr    Yda 
Ldr    Lda 
0      0 
Ndr   Nda]; 

a  =  inv(in)*an; 
b  =  inv(in)*bn; 
c  =  eye(size(a)); 
d  =  zeros(size(b)); 

t  =  0:05:15; 

%  calculate  sideslip  per  bank  angle 

f=[00-l*CL]'; 

k  =  [  CnB  Cndr  Cnda 

C1B  Cldr  Clda 

CyB  Cydr  Cyda]; 

perphi  =  inv(k)*f; 
betaperphi  =  perphi(l); 

PHI=10*pi/180; 
BETA  =  betaperphi*PHI; 
xin  =  [  BETA  0  Pffl  0  ]' 
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[x,y,t]  =  initial(a,b,c,d,xin,t); 

plotO.yC.lV-'fcgri^gtextCbeta') 

%title('Dutch  Roll  Response  for  Sideslip  I.  C.s') 

hold  on 

plot(t,y(:,3V-,);gtext('phi,) 

hold  off 

xlabel(Time  in  Seconds') 

ylabel('Amplitude') 
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4.        roll.m 

%  Solves  the  full  4x4  lateral-directional  for  roll  response 
%  (transfer  function  gain  and  phase)  due  to 
%  a  harmonic  aileron  input 

clear 
dderiv 

in  =  [Ufps      0         0        0 
0       1         0    -l*Ixz/Ixx 
0       0         10 
0  -l*Ixz/Izz    0        1] 

an  =  [YB     Yp    g*cos(thetanaut)     Yr-Ufps 


LB     Lp          0 
0      1            0 
NB     Np          0 

Lr 
0 
Nr] 

bn=[Ydr    Yda 
Ldr    Lda 

0      0 
Ndr   Nda] 

a  =  inv(in)*an 
b  =  inv(in)*bn 
c  =  eye(size(a)) 
d  =  zeros(size(b)) 

w  =  logspace(-l,2,150); 
[mag,phase,w]  =  bode(a,b,c,d,2,w); 

fori  =  l:length(w) 
forj  =  1:4 
if  phase(ij)  >  0 

phase(ij)  =  phase(ij)  -  360; 
end 
end 
end 

clg 

flgure(l) 

%  Plotting  roll  angle  gain 


77 


semilogx(  w,  mag( : ,  3  ));grid 
%title('Harmonic  Aileron  Input') 
xlabel('Frequency,  rad/sec');ylabel('Response  Gain') 
pause 

figure(2) 

%  Plotting  roll  angle  phase 

semilogx(  w,  phase( : ,  3  ));grid 

%title('Harmonic  Aileron  Input') 

xlabel('Frequency,  rad/sec'),ylabel('Phase  in  Degrees') 
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